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Abstract. We extend and correct a recently proposed maximum-likelihood halo-independent 
method to analyze unbinned direct dark matter detection data. Instead of the recoil energy 
as independent variable we use the minimum speed a dark matter particle must have to 
impart a given recoil energy to a nucleus. This has the advantage of allowing us to apply 
the method to any type of target composition and interaction, e.g. with general momentum 
and velocity dependence, and with elastic or inelastic scattering. We prove the method and 
provide a rigorous statistical interpretation of the results. As first applications, we find that 
for dark matter particles with elastic spin-independent interactions and neutron to proton 
coupling ratio fn/fp = —0.7, the WIMP interpretation of the signal observed by CDMS-II-Si 
is compatible with the constraints imposed by all other experiments with null results. We 
also find a similar compatibility for exothermic inelastic spin-independent interactions with 

fn/fp = - 0 . 8 . 
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1 Introduction 

We do not know what the dark matter (DM), the most abundant form of matter in the 
universe, consists of. Weakly interacting massive particles (WIMPs) are the most extensively 
studied DM particle candidates, not only because of their theoretical appeal but also because 
they could be detected in the near future. 

Direct searches attempt to measure the energy WIMPs might deposit when interacting 
within a detector. Three direct search experiments, DAMA [1], CoGeNT [2-6], and CDMS- 
II-Si [7] have potential DM signals, while all the other searches have produced only upper 
bounds on scattering rates and their annual modulation [8-18]. 

The conventional analysis of direct search data relies on a specific model of the DM halo 
of our galaxy, often chosen to be the standard halo model (SHM). A halo-independent data 
comparison method has been developed more recently [19-24]. The basic idea behind this 
method is that all the dependence of the scattering rate on the halo model, in any detector, 
resides in the same function which we call of the speed V m\u and the time t. 

Since is common to all direct search experiments, this function can be mea¬ 

sured by all experiments, and the compatibility of the different measurements can be studied. 
The speed Umin is the minimum speed necessary for the incoming interacting DM particle 
to impart a recoil energy En to a nucleus in each detector. Conversely, given an incom¬ 
ing WIMP speed v = Uminj is the extremum recoil energy (maximum energy for elastic 
collisions, or either maximum or minimum for inelastic collisions) that the DM particle can 
impart to a nucleus. Notice that En and Umin are exchangeable variables only for a single 
nuclide. When a target consists of multiple nuclides, a choice must be made between the two, 
Er and Umin- Taking Er as independent variable (as is done in [19, 20, 22]) Umin depends on 
each target nuclide. In our approach, Umin and the observed energy E' are the independent 
variables. This allows us to incorporate any isotopic composition of the target by summing 
over target nuclide dependent £'ij(umin) for fixed observed E'. 
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In earlier implementations of the halo-independent method, only weighted averages over 
nmin intervals of the time average, ti^(vmm), and annual modulation amplitudes, ?7^(umin), of 
vivmin,t) have been obtained from putative DM signals in direct detection. These averages 
over Umin intervals are represented in plots by a set of crosses in the Umin — V plane, whose 
vertical and horizontal bars show the uncertainty in 7?°(umin) or and the Umin range 

where they are measured, respectively. Combined with upper limits, these crosses can be 
used to assess the compatibility of data sets from various experiments. However, making 
a statistically meaningful evaluation of the compatibility of the data in this manner is not 
possible. 

The compatibility of different data has been studied in [25] using the “parameter 
goodness-of-fit” test statistic [26]. The analysis is based on the likelihoods maximized with 
fj^ written as a sum of very large number of step functions, following a method presented in 
[27]. In this case, the level of compatibility is given by the p-value of the test statistic, which 
was calculated by Monte Carlo simulations in [25]. Another test statistic for comparing one 
data set with a positive result and another with a negative result has been defined in [28]. 

An alternative method to study the compatibility of a positive result with upper limits 
uses a band in Umin — space at a given confidence level [29], derived from unbinned data, 
with an extended likelihood [30]. In this case, as shown in [29] for single-nuclide detectors, 
the likelihood is maximized by a non-increasing piecewise constant fj^ function, because of 
the exponential prefactor in the extended likelihood. 

The proof presented in [29] relies on the assumption that the target is made of a single 
component. The main limitation of the approach of [29] relies on their use of the recoil 
energy Er as independent variable. Here we provide a derivation of the extended maximum 
likelihood halo-independent (EHI) analysis method using Umin as a variable, which applies 
to any type of WIMP interaction, including inelastic scattering, and any target composi¬ 
tion. We correct and extend the original proof of [29] using the formulation developed for 
the generalized halo-independent analysis in [24]. The proof for the realistic case of finite 
experimental energy resolution presented in [29] relies on the application to the likelihood 
functional maximization of the Karush-Kuhn-Tucker (KKT) conditions in [31, 32]. The 
KKT conditions in [31, 32], however, apply only to the minimization of functions with a 
finite number of variables subject to a finite number of inequalities, and they do not apply 
to functionals. Eqs. (A.3) to (A.6) of [29] are given without proof and without a reference. 
Moreover, Eq. (A.4) seems problematic for a g function (which in our paper we call g) that 
has discontinuities, as in the solutions found in [29]. In this case, Eq. (A.4) requires a Dirac 5 
function to be smaller or equal to zero, which is mathematically problematic. As we explain 
in Sec. 3, although the KKT conditions have been extended to functionals defined on specific 
kinds of function spaces and constraints, we did not find in the literature a proof that clearly 
applies to our problem. Thus, in Sec. 3 we present our own proof of the KKT conditions we 
use, Eqs. (3.22)-(3.25), which are clearly valid for discontinuous functions. 

As in [29], here we find that the best fit fj function is piecewise constant with a number 
of discontinuities at most equal to the number of observed events. In [29], this is a result 
found for g given as a function of the recoil energy, which can be easily translated to Umin 
space only for a single target nuclide. Besides, the proof in [29] applies only to resolution 
functions with certain properties. We instead prove the result for fj as a function of Umin for 
any target composition and general resolution functions. 

Besides these extensions, we make a correction to the method of [29] by providing a 
clear definition of the uncertainty band. In [29], the uncertainty band is defined in Eq. (2.16) 
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through a numerical Monte Carlo simulation. In Sec. 4, we explain our objections to this 
procedure. We define instead a pointwise confidence band with a new method (see Sec. 4.2) 
and provide a clear statistical interpretation for this band using Wilks’ theorem. The different 
definitions of the band, here and in [29], yield very different values of the parameter AL 
defined in both papers for the same confidence level. 

In Sec. 2, we review the formulation of the generalized halo-independent analysis, on 
which the following sections are based. In Sec. 3, we prove crucial properties of the extended 
likelihood for unbinned direct dark matter detection data. In Sec. 4, we develop the EHI 
analysis method, and discuss the statistical interpretation of the confidence band computed 
with this method. In Sec. 5, we apply the method to the CDMS-II-Si [7] data for WIMPs 
with elastic isospin-conserving and isospin-violating SI interactions [33-35], and exothermic 
inelastic isospin-violating SI interactions [36, 37], and compare the results with the upper 
limits imposed by other experiments. Finally, we give our conclusions in Sec. 6. 


2 Generalized halo-independent analysis method 


The differential recoil rate per unit detector mass, typically given in units of counts/day/kg/keV, 
for the scattering of WIMPs of mass m off a target nuclide T with mass tut is 


(IRt 

dEji 


P Ct 

7Tl Ul'J' J 


_ ‘^min 


d^u f{v,t)v^^{ER,v), 

{Er) 


( 2 . 1 ) 


where Ct is the mass fraction of nuclide T in the detector, Er is the nuclear recoil energy, 
p is the WIMP local energy density, f{v,t) is the WIMP velocity distribution in Earth’s 
frame, dar/dER is the WIMP-nucleus differential scattering cross section, and U min fEp) is 
the minimum WIMP speed needed to impart to the target nucleus a recoil energy Er. The 
revolution of Earth around the Sun introduces an annual modulation of f{v, t). In detectors 
with more than one nuclide in their target, the total differential recoil rate is 


dii s dRT 
<^Er ^ dER 


( 2 . 2 ) 


Allowing for the possibility of inelastic scattering of an incoming WIMP with mass m 
into another outgoing WIMP with mass m' = m + S, for iJ,T\S\/m? <C 1, Vmm{ER) is given by 


(Er) 


1 

^/imjER 


TTItEr 

-hh , 

Pt 


(2.3) 


where pT = (fn mT)/(m + rtiT) is the WIMP-nucleus reduced mass. By inverting this equa- 

T i 

tion one obtains the minimum and maximum recoil energies E^’ (v) that are kinematically 
allowed for a fixed DM speed v, 




p^v'^ 

2mT 


1± W1 - 


25 

PTV'^ 


(2.4) 


The minimum possible value of Uminj thus also of v, for the interaction to be kinematically 
allowed is vj = \f25jp^ for endothermic scattering, and vj = 0 for elastic and exothermic 
scattering. This speed value corresponds to the point of intersection of the two Ej^ branches. 
We define vs to be the smallest of the vJ values among all nuclides T in the detector. 
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Most experiments do not measure the recoil energy Er directly. They measure instead 
a proxy E' for it, such as the ionization or scintillation signals, subject to experimental uncer¬ 
tainties and fluctuations. These are represented by an energy response function Gt{Eji, E'), 
which is the probability distribution for an event with recoil energy Eji to be measured with 
energy E'. Including the experimental acceptance e{Eji,E') from various experimental cuts 
and efficiencies, which, in general, is a function of both the recoil energy Eji and the detected 
energy E', the differential event rate in the detected energy E' can be written as 

^ eiER,E')GTiEii,E')^. (2.5) 


By inserting (2.1) into (2.5) we can express the differential event rate in the detected 
energy E' as a double integral 


^ = — <Er,E')Gt{Er,E') 

dE' rriT Jo 

x/ f{v,t)v^-{ER,v), ( 2 . 6 ) 

Jv>v^iUEn) ^Er 


from where changing the integration order and extracting a reference parameter a^ef from 
the cross section we get 


dR 

d^ 


f ^3^ 

^ Jv>Vs 


f{v,t) dJi 
V dE' 


{E\v). 


Here the function dJi/dE' is 


dU 

cL^ 


{E',v) 


E 


dJiT 


{E\v), 


(2.7) 


( 2 . 8 ) 


where dJiTj^E' is defined as 


dJiT 

d^ 


{E\v) 


^ r"" ^'’\ERe{ER,E')GT{ER,E') — ^{ER,v) iiv>vj, 
^ mr JeI;-(v) cF^ei<^ER 

0 if n < vj. 


(2.9) 

In (2.7) and (2.9) we have written explicitly a parameter Uref extracted from the differ¬ 
ential cross section to represent the strength of the interaction. This is preferably, but not 
necessarily, the WIMP-proton scattering cross section. 

Here we consider only differential cross sections (and thus also dJl/dE' functions) which 
depend only on the speed v = |n|, and not on the direction of the initial WIMP velocity 
V. The cross section depends only on v if the incoming WIMPs and the target nuclei are 
unpolarized and the detector response is isotropic, as is most common. In this case, one can 
write the differential event rate in a simpler form as 


dR 

d^ 


^ reip 
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V cW 
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( 2 . 10 ) 
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where F{v,t) ^ f dr 2 „ We now define the function as 


Thus, 


and (2.10) becomes 
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(TreiP F{v,t) _ dfl{v,t) 

m V dv ' 


dR 

Je' 



dri{v,t) dR 
dv dE' 


iE',v). 


( 2 . 12 ) 


(2.13) 


Using that fj{oo,t) = 0 (see (2.11)) and dR/dE'{E' ,vs) = 0 (since e’^ (vs) = e'^~^{vs) and 
the integrand in (2.8) is a regular function), the integration by parts of (2.13) leads to 


dR 



dUmin f]{Vynm,t) 


dn 

d^ 


{E , Umin), 


(2.14) 


where we choose to call Umin the integration variable because it makes obvious the physical 
meaning of ry as a function of Uniinj and where we define the “differential response function” 
dTZ/dE' of the detector as 


dn 

dW 


{E , Umin) 


d 

dVjnin 


dn 

d^ 


{,E , Ujnin) 


(2.15) 


Notice that dn/dE' is a function of the target dependent recoil energies, E'^^{vniin), which 
are functions of the independent variable Umin- It is clear that, in (2.14), all the dependence 
on the halo model is in the iy function which is independent of the experimental apparatus, 
and thus is common to all direct detection experiments. Therefore, by mapping the rate data 
into ry, it is possible to compare the different experimental results without any assumption 
on the dark halo of our galaxy. 

With dR/dE' given in (2.14), the energy-integrated rate over an energy interval [E[,E2] 


is 





(2.16) 


Direct detection experiments can measure the time-average R^^, 


and the annual modu¬ 
lation amplitude R^^, of the rate, R[e[,e^] — F[e' e'] ^[E' E'] cos(27r(t — to)/yr), with 

phase to- The only source of time dependence in the rate is ry, thus 


R 


[E[,E',] 


'vs 

roo 


dn 

d^ 


lin fy“(t'min) / dE' 

Je[ 

roo 

I dUmin fj ('^min)^[_Ej,£;'] ('I'min) 
J Vs 


(2.17) 


where a = 0 or 1, and fj{v,t) ~ ry'^(u) -|- cos(27r(t — to)/yr)- 

If the energy-integrated response function n\E' E'l(^min) for a given energy interval 

[E[,E'2l 


’^[E(,E'](^min) 





( 2 . 18 ) 
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is a well-localized function in a Umin range, the measurements of ^,j and g,j can be 

used to infer the values of fj^ and over the Umin range in which the response function is 
non-zero. This is true for WIMPs whose differential cross section is inversely proportional to 
such as for the usual spin-independent (SI) and spin-dependent interactions. Otherwise 
R[e[,e' 2 \ need to be regularized (see [24] for the details). 

For SI interactions, the WIMP-nucleus differential cross section can be written in terms 
of the effective couplings of the WIMP with neutrons and protons, fn and fp, as 

+ iAr - Zr)(U/U)f (2A9) 

where ap is the WIMP-proton cross section, fip is the WIMP-proton reduced mass. At and 
Zt are the atomic and charge numbers of the nuclide T, respectively, and F^{Ee) is a nuclear 
form factor, for which we take the Helm form factor [38] normalized to F^{0) = 1. Using 
(2.19) in (2.8) we obtain 

r <En,E')OT[ER,E') 

d-fc' ^ JeE-{v) 

x[Zt + (At - ZT)(fn/fp)?F^{ER), ( 2 . 20 ) 

and from (2.15), we get the following differential response function: dTZ^^/dE' 
j-^SI _ 

^[Zt + (At - ZT){fn/fp)f 

e(E]^+{v^in),E')GT(E]^+{v^in),E')F^(El’^(v^,n)) 
e(E'^~ (Vmm), E')Gt{E'^~ (Vmm), E')Ft(E'^~ (Vmm)) ■ 

V=Un-iin 

( 2 . 21 ) 


du 


For elastic scattering this reduces to 

(u„,in,S') = y^2v^u.—e{ER,E')GTiERiv^in),E') 
ah' rriT 

X^[Zt + (At - ZT)(fn/fp)?E^(ER(Vrr.^)). (2.22) 

f^p 

3 Piecewise constant f)(r;min) resulting from the EHI method 

Most direct detection experiments measure energy-integrated rates and/or their annual mod¬ 
ulation amplitudes in given energy intervals. CDMS-II-Si gives instead the recoil energies of 
three candidate DM events. Most halo-independent analyses of the CDMS-II-Si candidate 
events have chosen a binning scheme, which is arbitrary and may lose some of the information 
in the data [22-24, 36, 37, 39, 40]. 
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Ref. [29] has introduced a halo-independent analysis method without binning. The 
method relies on the fact that the extended likelihood [30] yields piecewise constant functions 
as solutions of the likelihood maximization. The extended likelihood for unbinned data can 
be written as 


No 

C[fi{v^^)] = H MT 

a=l 


dRtot 

d^;' 


E'=E' 


(3.1) 


For simplicity we use fj here for the time-average component of the fj function (we call it fp 
in previous sections). Here Nq is the total number of observed events, each with energy E'^, 
with a = 1,..., No- NElf]] is the total number of expected events within the energy range 
-^max] detectable in the experiment, which we write as a functional of the function 

flivmin)- 

poo 

NE[n]=NBG + MT dUminr/(Umin)'^[E' . .E' j('l'min), (3.2) 

I L mm" ma.xj 

Jvg 

where Nbg is the expected number of background events 


Nbg = MT 



dE 


/ dR-QG 
dE' ■ 


(3.3) 


Here MT is the detector exposure, dRiot/dE' is the total predicted differential event rate 


dRtot _ di?BG dR 
dE' ~ dE' ^ d^ 

_dRBG^r. .d7^ 

— + j dUniin^('t’min) ('I’min)) 


(3.4) 


and dRBG/dE' is the differential rate of the background events. Writing the rate in this form 
allows to take into account a non-trivial target composition (not included in [29]), through 
the differential response function dTZ/dE', defined in (2.8) and (2.15), or in (2.21) for SI 
interactions. 

Without fixing the halo model, the likelihood function in (3.1) is actually a functional 
of the fj function. If there is no uncertainty in the measurement of recoil energies, for a 
single target nuclide it was proven in [29] that the likelihood is maximized by a piecewise 
constant fj function with the number of steps equal to or smaller than the number of observed 
events, Nq- The proof for the realistic case with a finite energy resolution presented in [29] 
applies only to resolution functions with certain properties such as having a single local 
maximum, and relies on the application to the likelihood functional maximization of the 
Karush-Kuhn-Tucker (KKT) conditions in [31, 32]. The KKT conditions in [31, 32] apply to 
the minimization of functions with a finite number of variables subject to a finite number of 
inequalities. The proofs in [31, 32] do not apply to functionals. The likelihood C[fj] in (3.1) 
is instead a functional of fi{vmin) subject to an infinite number of inequalities, one for each 
value of Umin- The inequality given in Eq. (A.4) of [29], dfj/dEB in our notation, is actually 
an infinite set of inequalities, one for each Eb- 

The KKT conditions have been extended to functionals defined on specific kinds of 
function spaces and constraints. Banach spaces have been considered extensively (see e.g. 
[41]). However, the functions we are looking for, i.e. step functions, do not have derivatives 
everywhere unless interpreted as distributions, and the spaces of distributions defined on non¬ 
compact intervals like [0, oo) are not Banach (under the usual weak-* topology). More general 
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spaces, i.e. locally convex topological vector spaces, have been considered by Dubovitskii 
and Milyutin [42], As explained in the book by R. B. Holmes [43] (see pages 51 to 53), 
the Dubovitskii and Milyutin theory applies to constrained sets of functions that have non¬ 
empty interiors. However, the set of non-increasing functions such as fj has empty interior. 
A function iy G S' is in the interior of a set S if there is a neighborhood around it that belongs 
to the set, but for non-increasing functions there is no such neighborhood. This is because a 
non-increasing function always has a non-monotonic function arbitrarily close to it. 

Since we did not find in the literature any proof that clearly applies to our problems, 
in the following we present our own proof by first discretizing the variable Umin into a finite 
set of values so that the KKT conditions are applicable, and then taking the continuum 
limit at the end. Our proof is heuristic only in that it does not address the convergence of 
the limit. 

For convenience, let us define a different functional of ry, (—2 times the log-likelihood), 

L[fj\ = —2 In £[?)]. (3.5) 

With this definition, finding the ry function that maximizes the extended likelihood is equiv¬ 
alent to finding the function that minimizes L. To simplify the problem, we discretize the 
'I’min space into a set of A' -|- 1 positive variables = vs + i x Av with i = 0 , 1 ,..., K , 
where Av = — vs)lK with a large enough constant value. At the end, we will 

take AT —)> oo while keeping constant. 

With a A'-dimensional vector fj = (iyo, fji,..., fjK-i), we can define a piecewise constant 
function fi(vrain',fj) given by 

iy(Umin; fj) = fji if <in < ^^min < (3.6) 

Notice that there is no loss of generality of the iy(umin) considered, since any physically 
meaningful function is the limit of a sequence of piecewise constant functions as the number 
of steps tends to infinity. The corresponding L functional becomes a function of the vector 
V, 

/^(ry) = A[ry(umin;^)]- (3.7) 

With this discretization we can formalize the minimization of the functional A as a limit of 
the function minimization of fjs, and by doing so we can safely apply the KKT conditions. 

The KKT conditions for minimizing the function fiif]) under the constraints fji > fji+i 
on its variables (i.e. requiring the piecewise constant function iy(umin; fj) to be non-increasing) 
are the minimization conditions for the function 


K-l 

fiiv, ^ ^ fLiv) + qi{fji+i - fji) (3.8) 

1=0 

with respect to the variables fj and q = {qo,... ,qK-i) considered as unconstrained, and 
the supplementary conditions qi > 0 and (yj(7yj+i — fji) = 0. Written explicitly, the KKT 
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conditions are: 


(KKT I)a 

(KKT I)b 

(KKT II) 
(KKT III) 
(KKT IV) 


dfL 

dfji 

dh 

dfjo 


+ qi-i - qi = Q, ioi I < i < K 


- 90 = 0 , 


qi > 0 , 


1 , 


fji+i - i)i < 0, and 

qi{fii+i - fji) = 0, or equivalently, 

gj(r)j_|_i — fji)/Av = 0 (no summation imposed). 


(3.9) 

(3.10) 

(3.11) 

(3.12) 

(3.13) 


Choosing r) to be a unit vector fji along the ith component fji, the first term on the 
left-hand side of (3.9) and (3.10) can be written as 


_d_ 

dfii 


^ 9 -> 

fiiv) = Vi • -T^fLiv) = lim 
drj €^o 


fLiv + efji) - fiiv) 
e 


Using (3.14), we now have 


— f (') = lim ^ ~ v)] 

dfji e 


(3.14) 


(3.15) 


Using L[fj{vmm,V + ^Vi)] = L[fj{vai\n', v) + ef7(^'min; f?i)]) (3.15) Can be written in terms of the 
functional derivative of the L functional, 


A 

dvi 


fiiv) = liin 

e->0 


L[viVmm, V) + e^(^min; Vi)] “ L[v{v min; v)] 


e 

6L 


/ dVmin v{Vniin‘,Vi)- . 

lo OviVmin) 


(3.16) 

(3.17) 


From (3.6) one can easily see that the function fj{vmin',Vi) ™ (3-17) has a rectangular shape 
with value 1 between and and zero everywhere else. The summation over i from 
z = 0 to j < K' — 1 of the left-hand side of (3.9) and (3.10) is thus 


dfL 




dfL 


(3.18) 


2 = 0 

j roo 


5L 


A/ dvmin v{vmm,Vi)Y^i -T “ 9i- (3.19) 

t:/ Jo SviVmin) 


i=0 


Note that in the integrand (3.19) Y^i=oV{vram,Vi) = ^(An “ Vrain)0{vmin - vs), thus 


/•C 

5i 


dUmin V{v min ',Vi) 


6L 


SfjiVmin) 


9i = 


L 


min §Ij 

Sfiiv ■ ) 


(3.20) 


Using an interpolation function (?(umin) satisfying (/(uU^) = qj, we finally conclude, using 
(3.20), that the KKT I conditions (3.9) and (3.10) imply 


L 


H At 
^ ^V{v] 


- 9«in) = 0- 


(3.21) 
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By taking a large enough K with fixed, we can find an integer j such that is 

arbitrarily close to a given Umin value, if vs < Umin < Therefore, in the limit K ^ oo, 

and thus An —>• 0, we can write the conditions (3.21) and (3.11) to (3.13) for continuous Umin 
and fj variables: 

(I) qivmin) = dn (3.22) 

Jvs ^V{v) 


(II) 

Q'(^min) 

>0, 

(3.23) 

(III) 

Ve > 0, 

+ e) < ??(nmin), and 

(3.24) 

(IV) 

Q'(^min) 

i?('Cmin + e) - i?('(^min) „ 

iim -= 0. 

e ^+0 e 

(3.25) 


Note that although we write the conditions in terms of continuous variables, they should 
always be understood as a limit of the conditions for discrete variables. 

Two direct consequences of (IV) (3.25) are: i) ?)(nniin) can be discontinuous only at 
the points where gfu m^n ) vanishes, and ii) ?)(nmin) is constant in an open interval where 
gfu min ) / 0. If there is an open interval where g(nmin) is zero, within the interval, (IV) 
is trivially satisfied. Therefore, fj{vmin) is a piecewise constant function with discontinuity 
points where glu m^n ) = 0. Let us examine the possible zeros of the g(nmin) function. 

Using (3.1) and (3.5) in (I) (3.22), we get 


r'i'min 

qivmin) = 2 / dn 

Jvs 




6fi{v) 


I vs 


5fi{v) 


dE' 


(3.26) 


E'=E' 


In (3.2) Ne is given in terms of TZ, given in turn in (2.18), where dlZ/dE' is in (2.15). Using 
these equations, (3.26) becomes 


d.7Y 

g(n^in) = 2MT / dE' —{E'^v^m) 
Je' . dE' 


No r 

■ 2 E 


a=l 


dv 


6 dRtof 

6fi{v) dE' 


/ 

d-Rtot 

E'=E'al 

dE' 


(.3.27) 


J E'=E' 


We define: 


f^max d77 

avmin) = MT dE' —(i^', U^in). 


(3.28) 


Using (3.4) and (2.15), we can write 

6 dRiot 


dn 


'vs 


6fi{v) dE' 


, /d7^^ 

dn ( 

vs 


and 


— (-E^ ) ^mm) 


7a [i?] = 


)LtI 

= Ha{Vmm) 


E'=E'„ 


E’=E'„ 


dR, 


tot 


dE' 


(3.29) 


(3.30) 


E'=E'„ 
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Figure 1. i?o(^’min) (left panel) and ^{vmin) (right panel) for elastic isospin-conserving SI interactions 
and m = 9 GeV, for the three events of CDMS-II-Si. 


Replacing (3.28) to (3.30) into (3.27), we obtain 


^(^^min) 


No 

= 2C{Vrain) - 2Y, 

a=l 


Ha{Vmin) 

laiv] 


(3.31) 


In this equation, the only fj dependence is in 7a[i?]- The functions ^(umin) and Ha{vmm) do 
not depend on fj. 

Fig. 1 shows the functions F7a(umin) and £(u min ) for the three candidate events of CDMS- 
II-Si assuming an SI cross section with fn/fp = 1 and m = 9 GeV. In order to explain the form 
of these functions, let us first consider a simple situation where the target material consists 
of a single nuclide, or multiple isotopes of the same element, as in CDMS-II-Si. In this case, 
the integrands of the different terms dT-Lx/dE' in (2.9) contributing to dT-L/dE' in (2.8) are 
similarly localized in Epi for all nuclides T (for a fixed E'). Notice that these integrands 
are independent of v if dax/dEn is proportional to v~‘^. In this case, the v dependence of 
dTix/dE' is only in the integration range [E^~(v), E'^'^(v)]. If so, as v increases, this range 
covers more of the region in which the integrand is non-zero. Thus, dH/dE' grows with 
u in a certain range. When v is large enough for the integration in (2.9) to cover all the 
region in which the integrand is non-zero, dUx/dE' becomes constant, and so does dH/dE'. 
This explains the step-like functional form of Ha{v-[a\n) given in (3.29), which is dl-L/dE' with 
E' = Ea, as can be seen in the left panel of Fig. 1. 

Looking at (2.9) which defines dTix/dE' for each nuclide T, we see that the only de¬ 
pendence on E' of the integrand is in e{Epi, E')Gt{Er, E'). To compute ^(umin), we need 
thus a double integration, first in E^ to obtain dT-Lx/dE', and then in E', after summing all 
dUx/dE' contributing to dH/dE'. If we exchange the order of integration, performing the 
E' integration first, we see that as Er increases, for Er very small the integrand eGx will 
be zero within the E' integration range. Then, the non-zero portion of eGx within the E' 
integration range will increase, then be entirely contained, and then decrease and become 
zero again. Thus, the resulting integrand in Er will be slowly changing in the Er range in 
which it is non-zero. As Umin increases, the integration range in Er encompasses more of the 
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Figure 2. Same as Fig. 1 but for a fictitious detector with target material composed of equal mass 
fractions of Si and Ge (see the text), and showing in addition Haivmin)/'Ya (blue line, right panel) 


slowly varying integrand, resulting in a smoothly increasing function ^(umin), as shown in 
the right panel of Fig. 1. Once Umin becomes large enough for the integration range in Eji to 
cover all the non-zero part of the integrand in ^(umin), this function becomes constant (see 
Fig. 1). 

Let us return to study the discontinuity points of the best-fit fj function which happen 
at the zeros of q(u min ) given in (3.31). For elastic {6 = 0) or endothermic (5 > 0) scattering, 
there is a region at small Umin values where both Ha and ^ vanish. Looking at (2.9), when 
E^{v) is below the experimental threshold, the integrand, in particular the acceptance e, 
is zero, thus dH/dE' = 0. In this Umin region the condition (j'(umin) = 0 is trivially satisfied, 
and the shape of the best-fit fj function is undetermined. 

Changes in produce changes in 7a[?)]. For values of 7 a which make the second 

term of the right hand side of (3.30) large enough to reach the first term 2^(umin) from below, 
Q'(^^min) (see (3.31)) has non-trivial zeros where ^ and Ha are non-zero. The non-negativity 
of q{v) > 0 means that q{v) = 0 only when the monotonically increasing function ^ touches 
the step-like Ha/’ja function from above. Since Hafja has Nq steps, this can happen 
only at a number of Umin values smaller than or equal to Nq- Examples of these functions 
J2a Ha/lai and q will be shown below in Figs. 5, 7 and 8. 

To guess the generic shape of the f(u min ) and Hn{v m\n ) functions for differential cross 
sections whose WIMP speed v dependence is different from oc let us assume the differ¬ 
ential cross section for a given Er behaves as for large values of v. One such example 

is that of WIMPs interacting with nuclei through a magnetic dipole moment, where n = 2 
at large v (see (3.9) of [39]). In this case, from (2.9) one can easily see that the shapes 
of the functions Ha{v)/v'^ and ^{v)/v^ should be similar to those of Ha{v) and ^{v) for a 
differential cross section proportional to v~‘^. Therefore, the argument given above can be 
used for q{v)/v'^, whose zeros are the same as those of q{v), leading to the same conclusions. 

When the target consists of several elements, each Ha has multiple step-like features, 
one for each element. This is illustrated in Fig. 2 for a fictitious CDMS-II-like detector 
composed of equal mass fractions of Si and Ge. We see in the left panel of Fig. 2 that for 
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each of the three elements there are two step-like features in Ha- One may naively expect 
that because in this case there are 2No step-like features in Yla ^a/^a, the number of zeros 
of the function g(umin) would equally double. However, this is not the case. Because ^ and 
Ha are independent of ry, by changing fj and thus ja in general one can make at most one of 
the two steps per observed event in touch the function ^(umin) from below. Thus 

the number of zeros of g(umin) is still at most Nq- This can be seen in the right panel of 
Fig. 2. 

In summary, in this section we proved that the fj function maximizing the extended 
likelihood is a piecewise constant function with a number of steps smaller than or equal to 
the number Nq of observed events. 


4 EHI analysis in the Umin-space 


In this section we show how to find the solution to the maximization of the extended likelihood 
in the EHI method, in the Umin^h space. As shown in the previous section, the best-fit 
function, which we call ?7BF(r'min) from now on, is a piecewise constant function with at most 
Nq steps (note that in the statistics literature the subscript “ML” for maximum likelihood 
is usually used instead of “BF”). We will also find a statistically meaningful confidence band 
around fjnviv m^r. ). which we will define as a pointwise confidence band. 


4.1 Finding the best-fit function iyBF(r^min) 

The properties of the fj function maximizing the extended likelihood we have proven in the 
previous section can be utilized to find t^bf- We can define a function of 2No variables, 

V = {vi,V 2 , ■ ■ ■, vno) fj = (ryi, ? 72 , • • •, iyATo), specifying the positions and heights of the 
No steps, as a restriction of the functional L[fj]: 

= L[iy^^‘='^(umin;v,h)]- (4.1) 

The piecewise constant function fj^^o) jg defined as 


fj^^°\vram;V,V) 


fja if Va-l < Umin < Va, 

0 if VNo < Vrain, 


(4.2) 


where a = 1,..., Nq- Here we assume Umin and Ua’s are all larger than vs, and the constraints 
(3.24) fja < fjh for a > b are satished. Since the function fj cannot change after the last step 
and it must reach zero for large Uminj h must be zero for Umin > vnq- We do not specify the 
value of fj^^o) below the minimum vs since the event rate is independent of it. Notice that 
(4.2) requires the definition of vq. We define vq = vs for convenience. 

From these definitions and the theorem we have proven, we can easily obtain ryBF and 
Lmiii) the minimum value of the functional L[fj], by finding Fbf and ryBF that minimize 
so that 

iyBF('ymin) = (^'min; f^BF, ^Bf) (4.3) 

and 


Amin — fr[hBF(^^min)] — ^[iy^ (^minj ^BFj i?BF)] • 


(4.4) 


From the definition (3.1) of the extended likelihood function, we can write in a 

simple form as 


ANo) _ 
JL — 


j-Va 

2Nbg + ‘^MT E Va / 

a=l 


dUmin^fE' . .E' l('f^min) 

Lamin’ maxJ ' ' 


No 

-2^ In 
i=l 


MT^fla r du„,in^(^^min) +MT 

^ dE' HK' 


,dRBG , 
dE' 


(4.5) 


J E'=Ei 


with A^bg given in (3.3). Dehning the J\fa and Aiai functions of Va as 

rVa 

■^a{v) — MT / dViEmE'[E' ■ 

/ L min' maxj 

JVa-\ 

ATZ 

Adaj(u) = MT I dUjiiin , jp/ (l^min)) 

Jva-i E'=E'^ 

and the fixed constants bi 


we can write (4.5) as 


ANo) 

JL 


h = MT 


No 


di?BG 


dE' 


(4.6) 

(4.7) 

(4.8) 


E'=E', 


No 


21VbG + 2 ^ fjaJ^a — 2 ^ In 


a=l 


2 = 1 


No 


^ ^ Va-^d-ai 4“ bi 


.a=l 


(4.9) 


The minimization of the function of 2No parameters ui,... ,vnq, fji,... ,fiNQ, 

subject to the constraints 


vi > V 5 , (4.10) 

Vb —Va > 0 and fja — fjb A 0 for a < b, (4.11) 

can be done numerically using a global minimization algorithm. In the implementation, we 
express in terms of Inija and use Inija instead of fja as variables, since ija span many or¬ 

ders of magnitude. This also accounts for the fja > 0 constraints, leaving only the constraints 
in (4.10) and (4.11) to be enforced in the minimization. Note that in general minimization 
algorithms may attempt to evaluate the function in regions where the constraints are not 
satished, and in these regions the function is not well dehned, thus a fictitious function 

must be used that grows smoothly with the absolute value of the unsatisfied constraints in 
(4.10) and (4.11). 


4.2 Finding the confidence band 

In order to compare the 7 )bf we obtained with the upper limits imposed by other experiments, 
we need a way to represent the uncertainty in our determination of ?7 bf- This can be achieved 
by finding a region in the V m\n -f} space satisfying a certain statistical criterion, analogous to 
the confidence interval in the usual analysis with a fixed halo model. The region in the Umin“f? 
space which is densely filled by the family of all possible ?)(umin) curves satisfying 

AL[7)] = L[7)]-L^i„< AL*, (4.12) 
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with given AL*, is a natural candidate to examine. The condition in (4.12) defines a two- 
sided interval around t^bf for each r;min value, and the collection of those intervals forms a 
pointwise confidence band in Vrain~V space. From now on we will call it simply “the confidence 
band”. 

Conceptually, computing the confidence band is a straightforward procedure, but in 
practice, finding all the fj functions satisfying (4.12) and constructing the band from them 
is not possible. If the same band can be formed by a much smaller subset of them, and 
this subset is much easier to find than the whole set, the construction of the band would be 
practical. 

As a possible subset, let us consider the set of fj functions which minimize L[fj] subject 
to the constraint 

= (4.13) 

Let us define ,fi*) to be the minimum of the L[fj] subject to the constraint (4.13), 

and 

- ^min. (4.14) 

If ?)*) is larger than a chosen AL*, it simply means that the point {v*,fj*) lies 

outside of the confidence band. If it were inside the band, there should be at least one fj 
function passing through the (v*,fj*) point, for which AL[r)] < AL*, in contradiction with the 
fact that AL^-^^{v*,fi*) > AL*. On the other hand, if AL (j,-„ (n*, i)*) < AL*, the confidence 
band should cover the point {v*,fj*) by definition. Therefore, by finding the range of fj* 
values which satisfy AL^^^{v*,fj*) < AL* for each v* value, we can construct the band. 

The remaining problem is how to find an easy way of computing fj*) (and there¬ 

fore AL(j^-^(n*,?)*)). We will now prove that the fj function minimizing L[fj] subject to the 
constraint (4.13) should be a piecewise constant function with at most Nq + I discontinuities. 

Let us rewrite the KKT conditions in (3.9)-(3.13) but now with an additional equality 
constraint 

fjk = fj*, (4.15) 

where the index k is chosen to satisfy < v* < , so that can be arbitrarily close 

to V* for large enough K values. The additional constraint leads to the necessity of adding 
the term p*{fjk — fj*) to the function in (3.8) introducing a Lagrange multiplier p*, so we 
define another function f'l{fj,q,p*) as 

fLiv,Q,p*) = fLiv,^+p*im-fj*) (4.16) 

K-l 

= L[fj{Vmm,V)] + +P*iVk - V*), (4.17) 

1=0 

and use it to derive new KKT conditions. 

The new KKT conditions consist of the unconstrained minimization conditions of the 
function fi{fj,q,p*) with respect to the parameters fj, q and p*, plus the complementary 
conditions, which are the same as before. Therefore, besides the constraint (4.15), the changes 
only appear in the (KKT I)a and (KKT I)b conditions, where the function was present. 
The new conditions are 

-h qi-i - qi +p*5ki = 0, for 1 < z < iL - 1, and (4.18) 

dpi 

dfi 
dfjo 


(KKT I)'a 
(KKT I)'b 


qo + p*Sko = 0 


(4.19) 




with additional terms p*6ki and p*Sko, respectively, and the constraint (4.15). Following 
similar steps as those in Sec. 3 from (3.14) to (3.21), the summation of (4.18) and (4.19) over 
i from 0 to j now becomes 

/ du - u^) = 0. (4.20) 

In the limit of iF —)• oo, the first condition for the fj functions minimizing L[fj] subject to the 
constraint (4.13) becomes 

(I)' q{vrnm) = / dv +P*9{Vnim- V*), (4.21) 

Jvs O'qyv) 

while the conditions (II), (III) and (IV) are the same as in (3.23)-(3.25). 

Using the definition of L[ff\ in (3.5) and (3.1) in the condition (4.21), we can write the 
function qiv m\r, ^ as 


° R ( 1 

qiVrain) = 2e(^min) -2V +/0(u„,in - U*), (4.22) 

^ 7a W 

with ^(umin), Ha{vrain) and 7a[l?] defined in (3.28), (3.29) and (3.30), respectively. 

Again, the conditions in (4.21) and (3.25) tell that the fj function we find is piecewise 
constant with discontinuities only at the isolated zeros of q(vmm)- We already argued that 
£(v m\u ) can touch the function Yla=i from above at a number of points equal to or 

less than the number of observed events Nq- Since p*0{vmm — v*) introduces another step 
on the right hand side of (4.22), with the right p* value (?(umin) could have an additional 
zero. Thus the function minimizing L[fj] subject to the constraint (4.13) is piecewise 

constant with at most Nq + 1 discontinuities. 

Using a function fj of this type in (4.12) for each (v *we minimize L[fj] in (3.5) as in 
Sec. 4.1 to compute ,fj*) in (4.14). We define a function as in (4.1), 

parametrized by F = (ui, U 2 ,... ,Vi = v*,..., vno+i) and fj = {fji,fj 2 ,... ,fji = fj*,..., fjNo+i)- 
The minimization of can again be done numerically using a global minimization 

algorithm, subject to the same constraints as in (4.10)-(4.11), where in addition we keep 
{vi,fji) fixed at {v*,fj*). As before, in our implementation of the algorithm we write 
in terms of Inija instead of fja- We repeat the minimization procedure for all indices i = 
1... {No + 1) corresponding to the position of the {v*,fj*) step in fj, and select the solution 
that gives the overall minimum of 

4.3 Statistical interpretation of the confidence band 

From the procedure described above we can get both the best-fit fj function, i)BF('i’min), 
and the confidence band. For a quantitative assessment of the compatibility with other 
experimental data, we need to know the statistical meaning of a particular choice for AL*. 
One may be tempted to interpret AL as —2 times the logarithm of the likelihood ratio with 
2No parameters, since we parametrized the fj function with 2No parameters (plus v* and 
fj* which are fixed each time) to obtain the confidence band. However, this is not the proper 
interpretation. Note that the defining properties of the best-fit ryBF and the band do not rely 
on how we compute them. 
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Let us return to the definition of iy) and use again the discretization procedure 

introduced to derive the KKT conditions in Sec. 3. With a discretization of Umin we can define 
a likelihood function 

, VK-l) = A[77(Umin; V)] (4-23) 

with defined in (3.6). With this discretization, AL '^-^ {v*,fj) defined in (4.14) is 

replaced by a collection of functions each having v* in the /c-interval < v* < 

^mtn^ so that f]k = rf, 


Ai:;t(^*) = -21n 



. = V\Vk+l,---,VK-l) 


(4.24) 


Here the r)j values maximize the function C subject to the constraint i)fc = fj*, while % 
maximize the function £ without the constraint. Thus AL'^^ is —2 In of the profile likelihood 
ratio (see e.g. equation (38.53) of [44]) with only one parameter % = fj*. Notice that 
the continuous parameter v* becomes the discrete index k, and is no longer an additional 
parameter. According to Wilks’ theorem, the distribution of AL'^^^ approaches a chi-squared 
distribution with one degree of freedom, in the limit where the data sample is very large 
[44, 45] (and this is independent of the value of K). In short, this amounts to profiling 
the likelihood at fixed v* over the nuisance parameters fjo,..., • • ■, In this 

language, the fact that the likelihood ratio in (4.24) has one degree of freedom is proven 
mathematically in corollary 2 of [46] even for the case AT —)• oo. 

By taking large enough K, we can make and arbitrarily close to v*, and for 
each V*, AL'^^j^{fj*) approaches AL'^^^{v*,fj*). Therefore, the natural interpretation of the 
band is the collection of the confidence intervals in fj for each v^iin value, which defines a 
pointwise confidence band, based on a profile likelihood ratio with one degree of freedom. 
With this interpretation, we can now compare the confidence band with other limits or 
measurements in a statistically meaningful way. If any upper limit at some CL crosses the 
lower boundary of the band, at some other CL, it means that the two data, providing the 
limit and the band, are incompatible at their respective CLs. 

The Wilks theorem ensures the asymptotic behavior of the distribution of as 

the number of events becomes large, and the 3 observed number in CDMS-II-Si may not 
be a large enough number to ensure that AL follows the classical chi-squared distribution. 
Assuming that ALT^^ is chi-squared distributed, the choices of AL* = 1.0 and AL* = 2.7 
correspond to the confidence intervals of fj at the 68% and 90% CL, respectively, for each 
'I’min- The question on the convergence to the true confidence interval is also present in the 
analysis of the CDMS-II-Si data with a fixed halo model, if one uses the confidence interval 
estimator derived from the same likelihood function [22, 36]. 

In [29], AL* = 9.2 is used to compute the confidence band at the 90% CL, a value much 
larger than our choice, corresponding to the 90% CL limit for a chi-squared distribution with 
five degrees of freedom, resulting from a numerical Monte Carlo simulation. However, in 
the simulation described in [29], only fake data with three simulated events are selectively 
generated instead of allowing for any number of simulated events, as would be necessary to 
avoid generating a biased data set. Yet, allowing the number of simulated events to vary 
does not seem compatible with the AL definition in Eq. (2.16) of [29]. In this equation, V AL 
is defined as the radius of a hyper-ellipsoid in a 6-dimensional parameter space defined by 
the positions and heights of the three steps in the best-fit fjsF for a number of simulated 
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events Nq = 3. This leads to a chi-squared distribution for AL with 2No — 1 = 5 degrees 
of freedom (because there is one constraint). Allowing the number of simulated events Nq 
to change, the dimension of the fjBF parameter space is not fixed to 6, but would be 2No, 
leading to a number of degrees of freedom 2Nq — 1 that would change from simulated set to 
simulated set. 

5 Application to the CDMS-II-Si data 

In this section, we apply the EHI method to the three events observed by CDMS-II-Si in their 
signal region with recoil energies 8.2, 9.5, and 12.3 keV. We follow the procedure developed 
above. We use AL* = 1.0 and 2.7 for the 68% CL and 90% CL confidence bands, and compare 
the bands with the 90% CL upper limits from CDMSlite [16], SuperCDMS [18], LUX [17], 
XENONIOO [10] data, as well as the CDMS-II-Si data itself. The data analysis to obtain the 
upper limits in this paper is the same as in [36]. Recent analyses of the CDMS-II-Ge data 
[47, 48] use the same data set of [11], shown in [40] to provide weaker upper limits in the 
halo-independent analysis than SuperCDMS (and thus not included here). 

The data analysis described in this paper is implemented in the CoddsDM software [49] , 
an open-source Python program for comparing the data from direct detection experiments. 

5.1 Elastic SI scattering 

In this subsection we present the result of our analysis for elastic scattering with isospin- 
conserving fn/fp = 1 and with isospin-violating fn/fp = —0.7 (Xe-phobic) and fn/fp = —0.8 
(Ge-phobic) SI interactions [33-35]. They are shown in the left and right panels of Eig. 3 
and in Fig. 4, respectively, for a WIMP of mass m = 9 GeV. This value of the mass is 
within the 68% CL CDMS-II-Si regions obtained assuming the Standard Halo Model (SHM) 
in [40] and [36]. Figs. 3 and 4 show the best-fit fj^F (dark red line) and the 68% and 90% 
CL confidence bands derived from the CDMS-II-Si data shaded in darker and lighter red, 
respectively. Despite starting with three observed events, thus three steps in fj, the 7 )bf has 
only two steps, located at the zeros of the g(umin) function shown in Fig. 5. 

Fig. 5 shows the £(u min ) (red lines) and Xla-^a(^min)/7a (blue lines) functions in the 
left panel, and the g'(uniin) function given in (3.31) (right panel) for the best-fit fj^F of the 
CDMS-II-Si data for spin-independent elastic scattering with fn/fp = 1- The zeros of (?(umin)) 
located at the points where the functions in the left panel of Fig. 5 touch, are at 507 and 
580 km/s. These coincide with the locations of the steps of the best-fit j^bf plotted in the 
left panel of Fig. 3. The location of the steps is practically the same for other fn/fp values. 
The shapes of the £(^ 111111 )) Yla^a{vm\n)/la-, and q(u min ) functions are almost unchanged for 
a different choice of fn/fp values, up to a rigid rescaling along the vertical axis. The only 
changes expected in the positions of the zeros of g(umin) for different fn/fp values are due 
to the very small change in the relative strength of the WIMP interaction with different 
isotopes. Figs. 3 and 4 show the 90% CL CDMSlite (cyan), SuperCDMS (dark yellow), 
LUX (magenta), XENONIOO (blue) and CDMS-II-Si (red) upper limits, and the red crosses 
derived from the halo-independent analysis using binned data [40]. The crosses represent 
68% CL intervals of averaged fj and the corresponding Umin ranges for the CDMS-II-Si data 
with three equally spaced bins spanning the recoil energy range from 7 to 13 keV. Notice that 
the 68% CL crosses are similar in vertical extent to the 68% CL confidence band. Notice also 
that the 90% CL CDMS-II-Si limit follows closely the upper limit of the 90% CL confidence 
band. 
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Figure 3. 90% CL bounds from CDMSlite, SuperCDMS, LUX, XENONIOO and CDMS-II-Si, and 
the 68% CL and 90% CL confidence bands (see the text) from the three WIMP candidate events 
observed in CDMS-II-Si, for elastic isospin-conserving SI interaction {fn/fp = I, left panel) and for 
elastic Xe-phobic isospin-violating SI interaction {fn/fp = —0.7, right panel), for WIMP mass m = 9 
GeV. 



Figure 4. Same as Fig. 3, but for elastic Ge-phobic {fn/fp = —0.8) isospin-violating SI interaction. 


As one can see in the left panel of Fig. 3 for fn/fp = 1, the 68% CL confidence band 
is excluded in the Umin range from 370 to 560 km/s, by the combination of the 90% CL 
CDMSlite, SuperCDMS, and LUX upper limits. The lower boundary of the 90% CL con¬ 
fidence band is also cut at 450 km/s by the SuperCDMS 90% CL limit. Since there is no 
single continuous curve within the 90% CL confidence band which does not cross any 90% 
CL upper limit, we conclude that the potential signal and limits are incompatible for any 
halo model. 

On the other hand, in the right panel of Fig. 3 a significant portion of the 68% CL 
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Figure 5. ^{v^in) (red) and Y^a^i^o. {v-min)ha (blue) (left panel), and qivmin) = 2^{vniin) - 

2 J2a=i Haiy-cnin)ha (right panel) for SI elastic interaction with m = 9 GeV (see the text). 


confidence band remains below all the 90% CL upper limits. This shows that for SI interac¬ 
tions with fn/fp = —0.7 the CDMS-II-Si signal is consistent with the null results of all other 
experiments. 

The choice of fn/fp = —0.8 (Fig. 4) disfavors maximally the Ge limits (while fn/fp = 
—0.7 disfavors maximally Xe couplings instead). Thus, as expected, in Fig. 4 the SuperCDMS 
limit is weakened with respect to Fig. 3, but the LUX upper limits exclude almost completely 
both confidence bands. 

The dashed gray curves in Figs. 3 and 4 are the fj functions assuming the SHM for 
WIMP-proton cross sections Up = 10““^^ cm^ and ap = 10“^*^ cm^ in the left and right panels 
of Fig. 3, and ap = 10“^® cm^ in Fig. 4. For m = 9 GeV, these ap values are within the 68% 
and 90% CL CDMS-II-Si regions obtained assuming the SHM (in Fig. 1 of [40] and in Fig. 4 
of [36], respectively). In the analyses of [40] and [36] assuming the SHM, the m and ap choices 
for fn/fp = 1 and fn/fp = —0.8 interactions are shown to be rejected, while the choice for 
fn/fp = —0.7 interactions are allowed by all 90% upper limits. The same conclusions are 
evident in Figs. 3 and 4, where the dashed gray lines are above the upper limits for fn/fp = 1 
and fn/fp = —0.8 and below them for fn/fp = —0.7. 

5.2 Inelastic SI scattering 

In this subsection we present the results of the analysis for the exothermic Ge-phobic WIMP 
proposed in [36, 37] as an interpretation of the CDMS-II-Si data, shown in Fig. 6. This choice 
of fn/fp = —0.8 suppresses maximally the coupling to Ge. The limits due to Xe are weakened 
by the exothermic nature of the scattering, which disfavors heavier targets (such as Xe) with 
respect to lighter ones (such as Si) [36], leaving in principle Ge limits as the most important. 

Fig. 6 shows our results for a WIMP mass m = 3.5 GeV and mass split 5 = — 50 keV 
(left panel), and m = 1.3 GeV and 5 = —200 keV (right panel). These masses are shown 
in [36] to be within the CDMS-II-Si 90% and 68% CL regions when assuming the SHM, for 
ap = 10“^^ cm^ and 10“^^ cm^, respectively (see Figs. 5 and 6 of [36]). This is corroborated 
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Vmin [km/s] Vmin [km/s] 

Figure 6. Same as Fig. 3, but for Ge-phobic isospin-violating SI interaction {fn/fp = —0.8) with 
m = 3.5 GeV and S = —50 keV (left panel), and m = 1.3 GeV and (5 = —200 keV (right panel). 




Figure 7. Same as Fig. 5, but for SI exothermic inelastic interaction with m = 3.5 GeV and 6 = —50 
keV (see the text). 


by the present halo-independent analysis, since the corresponding fj functions assuming the 
SHM shown in Fig. 6 (dashed gray lines) escape all upper limits from experiments with null 
results. 

The best-ht ?7 bf functions for both Ge-phobic cases are shown in dark red in Fig. 6. 
They have two and one steps respectively in the left and right panels of Fig. 6, corresponding 
to the zeros of the g(umin) functions in the right panels of Figs. 7 and 8 (located at Umin 
values of 437 and 678 km/s in Fig. 7 and 792 km/s in Fig. 8). 

Figs. 7 and 8 show the functions ^ (red) and Yla=i Hn(v m\n )Hn (blue) in the left panels, 
and twice their difference, g(uniin)) in the right panels, for the two Ge-phobic cases in Fig. 6. 

In the previous analysis of Ge-phobic exothermic WIMP based on the SHM [36], the m 
and CTp parameters chosen in the current analysis are found to be compatible with the null 
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Figure 8. Same as Fig. 5, but for SI exothermic inelastic interaction with m = 1.3 GeV and 6 = —200 
keV (see the text). 


results of all other experiments. Consistently with this result, we find a large portion of the 
68% CL confidence band is below all the 90% CL upper limits imposed by all null results. 
Thus WIMP-nucleus scattering through Ge-phobic interaction can potentially explain the 
CDMS-II-Si data as a WIMP signal without any conflict with the null results of all other 
searches. 

6 Conclusions 

We have expanded and corrected a recently proposed extended maximum likelihood halo- 
independent (EHI) method to analyze unbinned direct dark matter detection data. Instead 
of the recoil energy as independent variable, we use Umiru the minimum speed a dark 
matter particle must have to impart a given recoil energy to a nucleus. An earlier version 
of the method, using Er as variable, was introduced in [29]. The use of Umin as variable 
allows to incorporate in the analysis any type of target composition and of WIMP-nucleus 
interaction, including elastic and inelastic collisions. This is not possible using Eji. The 
advantages of using Umin instead of Er in a halo-independent analysis was first pointed out 
in [21] and extensively used later on [23, 24, 36, 39, 40]. 

The EHI method uses unbinned direct dark matter detection data. The predicted 
differential rate as a function of the observed energy E' in all direct detection experiments 
can be written in terms of a common function r)(umin) (see (2.14)). The aim of the method 
is to find the fj function that provides the best fit for the unbinned data. We have proven 
rigorously that the best-fit fj function, ?7BF(^min)) is a piecewise constant function with a 
number of discontinuities smaller than or equal to the number of observed events Nq- We 
have also rigorously defined a two-sided pointwise confidence band with a clear statistical 
meaning, as a collection of confidence intervals in r) for every Umin value. We can assign a 
confidence level to the band and thus compare with upper limits given at particular confidence 
levels. This allows to quantitatively assess the compatibility of the unbinned data with upper 
limits due to null results. 
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Using this method, we analyzed the compatibility of the three candidate events found by 
CDMS-II-Si with the best available upper bounds, for spin-independent (SI) WIMP-nucleus 
interactions with different neutron to proton coupling ratio fn/fp values and either elastic 
or exothermic inelastic scattering. We found the best-fit t^bf function and 68% and 90% CL 
confidence bands. We chose values of the WIMP mass within the CDMS-II-Si regions in the 
m-ap plane that we had found in previous studies [36, 40] assuming the Standard Halo Model 
(SHM). Our results for fn/fp = 1, WIMP mass m = 9 GeV and elastic scattering are shown 
in the left panel of Fig. 3. The 90% CL CDMSlite, SuperCDMS and LUX limits derived as 
in [36] exclude the entire 90% CL band for this candidate. This case was also studied in [29], 
where the best-fit t/bf is very similar to ours, but the 90% CL band is much larger. In [29], 
only 90% CL limits derived from LUX and XENONIO data are presented. The LUX limit 
in [29] is similar to ours, but it does not exclude their much larger confidence band. 

The right panel of Fig. 3 shows our results for fn/fp = —0.7 (Xe-phobic) and m = 9 
GeV. We found that in this case a significant portion of the 68% and 90% CL confidence 
bands remains below all the 90% CL upper limits. Thus, a WIMP candidate with these 
characteristics provides an explanation for the three CDMS-II-Si events compatible with all 
present null results of other direct searches. This case was also studied in [29], where their 
best-fit 7 )bf function has the same number and position of steps as ours, but is an order of 
magnitude larger. We think this difference might be due, at least in part, to the inclusion of 
the isotopic composition of Si in our computation, which can not be done with the method 
used in [29]. The LUX limit presented in [29] for this case is similar to ours, but their 90% 
CL band is again much larger. 

The Ge-phobic fn/fp = —0.8 case, again for m = 9 GeV, is presented in Fig. 4. The 
90% GL confidence band is almost completely excluded by the 90% CL LUX limit. 

Our results for the Ge-phobic coupling and exothermic inelastic scattering are presented 
in Fig. 6, for two different values of the WIMP mass m and mass split 6: m = 3.5 GeV, 
S = —50 keV and m = 1.3 GeV, S = —200 keV. In these cases the 68% and 90% GL confidence 
bands are almost entirely below all the 90% GL limits. Thus, again we found compatibility 
between a dark matter interpretation of the CDMS-II-Si data and all null results. 

In all cases studied we included the crosses derived from the CDMS-II-Si data obtained 
with our previous halo-independent analysis using binned data [36, 40]. The crosses represent 
68% CL intervals of averaged rj and Vmin ranges corresponding to three equally-spaced bins 
spanning the recoil energy range from 7 to 13 keV. We found that the crosses are similar in 
vertical extent to the 68% CL confidence bands in all cases. This shows agreement between 
both types of halo-independent analyses, but the present method is much more powerful. 

We found remarkable that the 90% CL limit derived from the CDMS-II-Si data itself 
using the Maximum Gap method, as described in [36, 40] (and references therein) is almost 
identical to the 90% CL upper boundary of the 90% CL confidence band in all cases studied. 
Again, this indicates agreement between the two different analyses. 

SI elastic scattering was also studied in [25] and [28], where two different statistics 
were used to quantify the compatibility among different direct search data sets. In [25], for 
isospin-conserving SI interactions and WIMP mass 7 GeV, which is slightly smaller than 
our choice of 9 GeV, the parameter goodness-of-fit value derived from the global likelihood 
of the CDMS-II-Si, SuperCDMS and LUX data has a p-value of only 0.44%. This poor 
compatibility level is consistent with our results. For isospin-violating interactions, [25] used 
slightly different parameter sets, fn/fp = —0.71, m = 6.2 GeV, and fn/fp = —0.79, m = 6.3 
GeV, with corresponding p-values of 18.7% and 18.5%. Thus the compatibility is significantly 
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improved, which is also consistent with our results. In [28], a test statistic “pjoint” is proposed 
and calculated, said to be the upper bound on the joint probability of obtaining the outcomes 
of two potentially conflicting experiments. Only if the value of pjoint is small there is a clear 
interpretation of incompatibility, but a large pjoint value does not imply compatibility. For 
m = 9 GeV, [28] finds incompatibility between CDMS-II-Si and SuperCDMS for fn/fp = 1, 
but not for fn/fp = —0.7 or —0.8. In this respect, we agree. 

The use of a test statistic such as defined in [25] or [28] is complementary to our method 
of using a confidence band and upper limits in Vmm—V space to assess the compatibility among 
different data sets. 
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